Quasiclassical numerical method for mesoscopic superconductors: bound states in a 

circular d-wave island with a single vortex 
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We propose a numerical method for studying mesoscopic superconductors in terms of the qua- 
siclassical theory of superconductivity. Our method is free from ambiguity in initial values for a 
system that does not have a bulk solution, and allows one to uniquely determine selfconsistent 
solutions of the Eilenberger equation which are independent of initial conditions for quasiclassical 
trajectories. To demonstrate the efficiency of our method, we calculate the local density of states in 
a circular d-wave island containing a single vortex. We find that the so-called vortex shadow effect 
strongly depends on energy in a circular mesoscopic island. 

PACS numbers: 74.78. Na 74.20. Fg 74.25.Ha 



I. INTRODUCTION 

Recently developed experimental techniques have 
made it possible to fabricate mesoscopic superconductors 
and to observe their electronic structure by scanning tun- 
neling spectroscopy (STS).^!^ Due to finite-size effects, 
mesoscopic superconductors can exhibit properties that 
are significantly different from those of their analogous 
bulk materials. For example, vortex physics presents var- 
ious intriguing phenomena in a mesoscopic system whose 
size is of the order of the coherence length or the pene- 
tration depth. In particular, competition between the re- 
pulsive interaction among vortices, which tends towards 
formation of the Abrikosov vortex lattice, and quantum 
confinement effects results in a variety of vortex states 
that are unique to small systems. The signature of gi- 
ant vortices carrying multiple flux quanta'^ and that of 
"shell effects" of multiple vortices, where vortices arrange 
themselves conforming to the shape of the sample, have 
been detected in submicron Al disks i^i^ Depending on 
the size and shape of the system, a pair of vortex and 
antivortex can also be formed.^ STS can directly probe 
the local density of states (LDOS) in such novel vortex 
states. 

It is important to determine the phase of the super- 
conducting order parameter in unconventional supercon- 
ductors such as cuprates, heavy electron superconduc- 
tors, and iron-based materials. One of the important 
characteristics of unconventional superconductivity is the 
possibility of the existence of Andreev bound statesji"— 
When there is a sign change in the order parameter in 
momentum space as in d-wave superconductors, Andreev 
bound states can be formed if the quasiparticle feels the 
sign change by specular reflection at a surface. Andreev 
bound states can also exist where the order parameter 
changes its sign in real space, e.g., around a vortex. The 
formation of Andreev bound states is thus a key phe- 



nomenon that can reveal the fundamental nature of su- 
perconductivity. 

In unconventional superconductors, phase-sensitive 
phenomena can be manifest in systems where interfer- 
ence effects can occur between a vortex and a surface. 
In (ja;2_j,2-wave superconductors, the "vortex shadow" 
effect, which suppresses zero-energy density of states, 
has been observed near a vortex in front of a reflecting 
110-boundary.^" In chiral p-wave superconductors, low- 
energy Andreev bound states can be either suppressed or 
enhanced by a vortex, depending on its orientation with 
respect to the chirality of p-wave superconductivity»ii 
Such phase-sensitive phenomena are expected to appear 
in mesoscopic superconductors, where the effects of sur- 
faces can be dominant. 

The electronic structure of the vortex state has been 
studied in terms of microscopic mean-field theory, with 
which one can calculate the LDOS observable by STS. 
Such a microscopic calculation was first performed by 
Gygi and Schliiter,^^ who obtained the LDOS around a 
vortex by numerically solving the Bogoliubov-de Gennes 
(BdG) equations. With the use of the quasiclassical the- 
ory of superconductivity) "'^'^'"'^'^ Hayashi et alJ^ have repro- 
duced the LDOS in the vortex state observed in NbSe2 
by STS.i- . The electronic structure around a vortex in a 
d-wave superconductor has been calculated by Schopohl 
and Makifii using the Riccati parametrizationi^ of the 
Eilenberger equation in the quasiclassical theory. The 
Riccati formalism has also been developed by Ashida 
et al. in the context of a boundary problem for 
superconductor-normal-metal interfaces; 

Recent STS measurements have shown the direct evi- 
dence of giant vortices and multivortex configurations in 
nanoscale Pb islands. Due to limited resolution, how- 
ever, bound states corresponding to different numbers 
of flux quanta could not be identified in the measured 
LDOS in either experiments. Theoretically, the LDOS 
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in the giant vortex state in s-wave mesoscopic disks has 
been calculated by directly solving the BdG equationsi^ 
Solving the BdG equations, however, has high numeri- 
cal demand, while in contrast the Eilenberger equation is 
relatively easy to solve, especially by means of the Ric- 
cati parametrization. However, in order to integrate the 
Riccati equations one needs to know the initial values of 
the Riccati amplitudes, namely, the boundary values in 
the case of a finite-size system. For this reason, most 
of the systems studied so far have a bulk region, from 
which to integrate the Riccati equations with well-defined 
initial values. One exception is the work on the vortex 
lattice;^ which does not have a "bulk" region in the sense 
of having unique initial values. Miranovic et al. have 
calculated the electronic structure in the vortex lattice 
within the Riccati formalism, -5^. although they have not 
discussed why their method can be numerically stable for 
solving the Riccati equations. The vortex lattice systems 
have also been studied by solving the Eilenberger equa- 
tion directlyi^i^ To the best of our knowledge, there has 
seldom been work on rigorously solving the Eilenberger 
equation for a system which has no bulk solution such as 
finite-size superconductors. 



In this paper, we propose an efficient and stable nu- 
merical method for studying mesoscopic and nanoscale 
superconductivity within the framework of the quasiclas- 
sical theory. Our method releases one from the problem 
of determining the initial values in a system that does not 
have a bulk region, and makes it feasible to obtain initial- 
value-independent solutions to the Eilenberger equation. 
We have developed this method by generalizing the nu- 
merical technique of Ref. [2l| so that it is applicable to 
finite-size systems and suitable for systematic exploration 
of mesoscopic and nanoscale superconductors of various 
size and shape. We provide an explanation as to why 
this method is numerically stable in the Riccati formal- 
ism as well as the condition for stability. To illustrate our 
technique, we calculate the LDOS in a circular d-wave is- 
land sustaining a single vortex and find that the "vortex 
shadow" effect strongly depends on the quasiparticle en- 
ergy in such small islands. 



The paper is organized as follows. In Sec. II, we sum- 
marize the Riccati formulation of the quasiclassical the- 
ory of superconductivity and show how to obtain initial- 
value-independent solutions by numerical computation. 
In Sec. Ill, we introduce our model of a circular dx2-y2- 
wave island containing a single vortex, and present re- 
sults for this system and main conclusions in Sec. IV. 
and V, respectively. A general solution of a Riccati-type 
equation is presented in Appendix A. In Appendix B, we 
discuss the stability of the Riccati equations with the use 
of the analytical solutions of the bulk and in the vicinity 
of a single vortex. How to generate a path with specular 
refiections inside a circular disk is shown in Appendix C. 



II. FORMULATION 

A. Quasiclassical theory of superconductivity 

We introduce the quasiclassical Green function g de- 
fined by 



g{iujn,r,kF) 



9 f 

-f -9 



(1) 



which is a function of the Matsubara frequency ujn , the 
Fermi wave vector fcp, and the spatial coordinate r. The 
check A signifies the 2x2 matrix structure in the Nambu- 
Gor'kov particle-hole space. The Eilenberger equation is 
the equation of motion for g{iujm r, fcp), 

-ivpikp) ■ Vg = [iQnfs - A(r, fcp), g] , (2) 

supplemented by the normalization condition, 

(3) 

where iuJn = i^n + i'f • f ^ with A a vector potential and 
73 is the Pauli matrix. The A(r, fep) is given by 



A(r,fcF 



A(r,feF) 
-A*(r,A;F) 



(4) 



in the Nambu-Gor'kov space. Setting ia;„ — e-\-ir], where 
rj is real and positive, we obtain the retarded quasiclas- 
sical Green function. 



B. Riccati formalism 

We solve Eq. ^ by means of the Riccati parametriza- 
tion. 



-1 
1 + ab 



l-ab 
-2ib 



2ia 
(1 - ah) 



(5) 



Using these variables, the Eilenberger equation ([2|) re- 
duces to a set of two decoupled differential equations of 
the Riccati type. 



vp ■ Va = — 2a)„a — A*a^ 
t>F ■ V5 = +2ujnb + A6^ - 



(6) 
(7) 



Since these equations contain V only through tJp ■ V, 
they can be reduced to a one-dimensional problem on a 
straight line in the direction of the Fermi velocity vp, 

da. . u, o . 

= -2una - A*a^ + A , (8) 



VF 



ds 
db 



-2Lj.„b + Ab^ 



A* 



(9) 



The local density of states (LDOS) for an isotropic Fermi 
surface as a function of quasiparticle energy e (with re- 
spect to the Fermi level) is given by 



dfl]s 

in 



■Rc 



1 - ab 
1 + 06 



(10) 
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where i^(0) is the Fermi-surface density of states, dJlk is 
the solid angle, and is a smearing factor as due to im- 
purity scattering. We now demonstrate how to obtain 
initial-value-independent solutions to these Riccati equa- 
tions without the need of setting initial values. 



C. Initial-value-independent solution 

Let us consider the Riccati equation with complex fre- 
quency z: 



da . » T . 

vf— = 2iza - A*a^ -\- A . 

OS 



(11) 



If we can find a particular solution a = ap(s), we can 
obtain a general solution as (See, Appendix A) 

a{s) = ap(s) + 



/;^^(s')e-^(^'))e^(-)+u(so) ' 



with 



Ais) = - 



(12) 



(13) 



K{s) = A r A*(s')ap(s') - 2i—is - so) . (14) 

VF Js„ VF 



The u(so) satisfies an initial condition at s = sq: 

1 



a(so) ap(so) -|- 



u{so) ' 



If the condition 



lim K{s) — -\-oo, 



(15) 



(16) 



From Eq. ([14]), if J^^ Re [A*(s)aN(s)] is an increasing 

function of s in the upper half plane of 2, e^^*') increases 
with increasing s, since the second term in Eq. (jl4p is 
always a monotonically increasing function in the upper 
half plane of z. The length is characterized by the fcp- 
dependent coherence length ^(fc) = vpikp)/ A{kF)- In 
the region s — sq 3> ^(fcp), we obtain 



(21) 



since the second term in Eq. (|19p will approach zero 
within a numerical accuracy. Therefore, one can always 
obtain a numerically stable solution aN(s), which means 
that the solution does not depend on the initial value if 
far away enough from the initial point. These discussions 
are originated from the fact that the numerical solution 
aj^{s) can be regarded as a particular solution. We find 
that the relation Re (A*(s)aN(s)) > is usually satisfied 
within a wide range in various systems (See Appendix 
B). This argument is valid when integrating the Riccati 
equation for b in Eq. 

We are thus released from the problem of deciding on 
the initial value in a system which does not have a bulk 
solution. Our discussion above clearly shows the rea- 
son why one has to integrate Eq. ([8]) in the direction 
of increasing s and Eq. ^ in the opposite direction of 
decreasing s. In the upper half plane of z, the second 
term in Eq. (|14p increases monotonically with increasing 
s. On the other hand, one has to integrate Eq. ([8]) in 
the direction towards decreasing s when considering the 
lower half plane of z. Our method is appropriate for the 
vortex lattice and analogous to the method of Ref. [2l| 
of integrating the Riccati equations for obtaining stable 
solutions. 



III. NUMERICAL METHOD 



is satisfied in the upper half plane of z, the solution a(s) 
does not depend on u{so) in the limit of s 00: 



lim a{s) — ap(s). 



(17) 



Now, we assume that we have obtained a numerical 
solution aN(s) with the initial value at s = sq, 



(18) 



If one wants a solution a^(s) with another initial value 
a'o at s — Sq, one can obtain it by 



a'-^is) = aN(s) + 



/; A(s')e-^(^')) e^(^) +u(so) 



where 



u{so) 



Aq - ao 



(19) 



(20) 



We illustrate our method by calculation of the LDOS 
in a circular d-wave island with a single vortex. 



A. Model 

We consider a two-dimensional system of circular 
shape of radius Tc, which has a specular surface and a 
circular Fermi surface. The boundary condition can then 
be expressed asiS 

a{\r\ = rc, fcm) = a{\r\ = r^, feout) , (22) 
6(|r| = r„fci„) = 6(|r| = r„fcout). (23) 

Here k-m is connected with feout by specular reflection. 
We introduce a pairing potential of the form. 



A(r,fcp)- Ao/(r)d(fcp)e' 



(24) 



where r = r(cosQ;, sina) in polar coordinates. Here /(r) 
gives the spatial variation of the pairing potential with 
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/(O) = at the vortex center, and d^kp) is the gap 
anisotropy in momentum space. The direction of the 
quasiparticle motion is characterized by the angle 6 in 
momentum space. We consider a c?2,2_j^2-wave supercon- 
ductor: 

d(A;F) = COS26'. (25) 

Considering the strongly type-II limit, we neglect the 
vector potential iw„ — >■ By setting iujn e + ir/, 

we integrate the Riccati equations by the fourth-order 
Runge-Kutta method. For simplicity, we define a gap 
function without solving the selfconsistent gap equa- 
tion. Throughout this paper, we use the length unit 
^0 = i'f/Ao, smearing factor 77 = O.OIAq and the spatial 
variation /(r) = r/ y^r"^ + £,q- In order to obtain a{x, y, 9) 
and b{x, y, 9), we must integrate the Riccati equations on 
the paths with specular refiections inside the circle. 



B. Numerical recipe 

We now describe how to obtain the values a and b at 
(xa,yo) with the momentum direction 9. At first, we 
generate path I with specular refiections inside the circle 
for a as shown in Fig.[TJ by drawing the path from (xq, yo) 
to {xn,Vn) (See Appendix C). The (x„,?/„) (n > 1) in 
Fig. [1] is the point at the n-th specular reflection. We 
generate path II for h in the opposite direction. It should 
be noted that the length of two paths after n reflections 
should be much longer than the coherence length. This 
means that the smaller the system, the larger n should 
be. 




specular reflection 



FIG. 1. (Color online) Schematic figure of the specular reflec- 
tions. 

Next, we integrate the Riccati equations for a (b) from 
(xn.yn) to (xo,2/o) (from {x'„,y'„) to (xo,2/o)) and obtain 
the values a and b at {xo,yo) in the direction 9. With 
the use of Eq. ((TU)) we can obtain the LDOS of a circular 
d-wave island. One must make sure that the results do 
not depend on the initial value or the length of the path. 



IV. RESULTS 

We present the LDOS with and without a vortex in a 
circular d-wave island. 



A. Without a vortex 

First we show the LDOS for a system without a vortex. 
We have used 720 0-meshes in momentum space. In the 
LDOS presented in Fig. [21 the zero-energy bound states 
can clearly be seen at 110-boundaries. This confirms that 
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FIG. 2. (Color online) Local density of states in a circular 
dj.2_j^2-wave island without a vortex at energy (a) e = 0, 
(b) O.O5A0, and (c) O.IAq. The radius is rc = 5^o and the 
smearing factor rj = O.OIAq. 

our numerical procedure is appropriate for small d-wave 
systems. 

B. With a single vortex 

Next we consider a single vortex at the center of the 
circular island. We can see in Fig. [3] the "vortex shadow" 
effect which has been discussed by Graser et al. at a sur- 
face of a d-wave superconductor At zero energy, there 
is no surface bound state because of the vortex shadow 
effect as can be seen in Fig. ^a.). Surface bound states 
exist with nonzero energy and their pattern changes with 
increasing energy. A trajectory in the region where the 
LDOS intensity becomes stronger near the vortex can be 
regarded as a "ray of light" for the surface bound states 
as shown in Figs.[31[b) and[3l^c). We demonstrate in Fig.|3] 
that the vortex shadow effect becomes small with increas- 
ing rc- This Tc-dependence is similar to the dependence 
of the distance between the vortex and the boundary in 
Fig. 5 in Ref. M- 
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FIG. 3. (Color online) Local density of states in the circu- 
lar dj.2_y2-waye island with a vortex at the energy (a):e = 0, 
(b):0.05Ao and (c):0.1Ao. The radius is Vc = 5^o- The smear- 
ing factor is 77 = O.OIAq. 



CONCLUSION 



In conclusion, we have proposed an efficient numeri- 
cal method for studying superconductivity in mesoscopic 
and nanoscale systems in terms of the quasiclassical 
Eilenberger equations. We have described the numerical 
procedure by which a stable, initial-value-independent 
solution can be obtained, due to the fact that a numer- 
ical solution can be regarded as a particular solution of 
the Riccati equation. Our technique has been illustrated 
with the LDOS in a circular d-wave island with a single 
vortex. We find that the "vortex shadow" effect strongly 
depends on the quasiparticle energy in such small sys- 
tems. For the sake of illustration, we have assumed the 
spatial variation of the gap function. It is straightfor- 
ward, however, to incorporate our method into solving 
the gap equation selfconsistently, as well as to include a 
vector potential in selfconsistent calculation. 
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FIG. 4. (Color online) Local density of states in the cir- 
cular d^2_y2-wave island with a vortex at the boundary 
r = r-c(cos a, sin a) as a function of polar angle a and energy 
e, for radius (a) = 5^0, (b) lO^o, and (c) 20^0. 



Appendix A: General solution of the Riccati 
equation 

We consider the Riccati equation in the general form, 

(Al) 



ax 



If we can find a particular solution y = f{x), we can then 
obtain a general solution using the form y — f{x) + l/u. 
The differential equation for m is a linear equation, 

nil 

^^~i2A{x)f{x)+B{x))u~A{x). (A2) 
ax 

When we consider the initial value u(xo) at xq, the solu- 
tion can be expressed as 

u{x) = - (^j" da;'A(a;')e"^("')^ e^(") + u{xa), (A3) 
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with 

K{x) = - I dx' {2A{x')f{x') + B{x')) . (A4) 

We thus find the general solution as 

1 



y = fix) + 



J!^^ dx'A{x')e-K e^(^) + u{xn) 



(A5) 



Appendix B: Stability of the Riccati equations 

We discuss the stability of the Riccati equations with 
use of the analytical solutions. 

1. Bulk 

The solution for a homogeneous bulk system is 



-Lj„ + via! 



A* 



Therefore, Re(aA*) > when e < A. 



(Bl) 



2. Near a vortex 

Near a vortex, we can use the Kramer-Pesch approxi- 
mation (KPA);^"^ The KPA can be regarded as a per- 
turbation with respect to both the quasiparticle energy 
and the imaginary part of the pairing function in the 
Riccati formalism. Introducing the variables, 

(B2) 
(B3) 
(B4) 



a — 






be-' 


A = 


Ae'« 



the Riccati equations can be rewritten as 

d _ ^ ^ _ — 2a* a 

Dp— a = —zojna — a A -f A , 

OS 

d - 

VF^b = 2cj„fe + b^A - A*. 
OS 



(B5) 
(B6) 



The two-dimensional polar coordinates are denoted here 

as 



r = (s, y) = r(cos 6, sin 6'). 
In these coordinates, A reduces to 

s + iy 



A(r,feF) =/(r)Aod(fcF)- 



By means of KPA, we have 

d{r, fcp) ~ ao(A;F) -I- ai(r, fcp), 



(B7) 



(B8) 



(B9) 



where 

Oo(fcF) 

ai(r,fcF) 
Here, 



= -sgn [dikp)] 
= -2- 



VF 



(BIO) 

[aoikppn - ilmA(r')] e^^^'^'^ds'. 

(Bll) 



u{ 



r) = — ao(feF) / ReA{r')ds'. (B12) 



The condition Re(A*a) > is then expressed as 

2g«(r) 



ReA*a = D(s) 



VF 



{sio„Cis) - y^Eis)) 



with 



C'is) 
D{s) 
E{s) 



f{r)Ao\d{k)\ 
y^s^Ty^ 

D{s')e-''^'''Us'. 



>0, 
(B13) 

(B14) 
(B15) 
(B16) 



Since e"*^''-' is a localized function at s = and the ap- 
plicable range for the perturbation is |ao| > |ai|, the 
condition Re(A*a) > is satisfied in the region s < 
and ujn > 0. This means that one can obtain numer- 
ically stable solutions in a system containing a vortex. 
Furthermore, as K{s) is a function obtained by the inte- 
gration of Re(A*a), it is an increasing function close to 
and far away from a vortex, and can yield numerically 
stable solutions in a system with many vortices as long 
as the intervortex distances are sufficiently longer than 
the coherence length. 



Appendix C: Path with specular reflections inside a 

circle 



We illustrate how to generate a path with specular re- 
flections inside a circular disk, from a initial point (xq, yo) 
with the initial direction 9. The linear path that goes 
through the point [xq, yo) with the gradient a — tan^ is 
written a,s y = a{x — xq) + yo. We flnd the point of in- 
tersection of this path with the circular boundary, which 
is given by x'^ + y'^ — r'^. The solutions are 



x± 

y± 



a?XQ — aya ± D 
a{x± - Xo) + yo , 



with 



D = \Jrl+ a^rl - a'^x^ + 2axoyo - yo 



(CI) 
(C2) 

(C3) 
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Denoting {xc,yc) = {x-tU-) as the point of intersection, 
we can obtain the path as 



y ^ a{x - xq) + uq, {xc<x<xo). (C4) 

The angle of specular reflection 0' can be found by simple 
geometry: 



We obtain the next path with the direction Oi-i 
by adopting as the initial point on the z-th path 
{xi-i,yi^i) = {xcVc)- Then, we find the point of in- 
tersection: 



x± 



a^g;,,_i - ayj-i ± D 
l + a2 

y± = a{x± - x,_i) + , 



(C7) 
(C8) 



9' = e + 269 



a- 9 (a > 0) 
t: - 6 + a (a < 0) 



(C5) 
(C6) 



where (xc,yc) = r(cos a, sin a) in polar coordinates. We 
thus regard the angle 9' as the new momentum direction 



D 



\Jrl+ a?rl - a?x1_^ + 2aa;j_ii/i_i - y'j_^ , (C9) 



with a ~ i'ATi9i-i. One of the solutions is equal to Xi-i 
and the other solution should be adopted as the next 
intersection point. Then we obtain the z-th path: 



y = a{x ~ Xi^i) +yi^i, {xi^i < x < Xi) . 
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